*12345678901234567890123456789012345678901234567890123456789012345678901234567890
capture log close
clear all
set more off

*	************************************************************************
* 	File-Name: 	Aklin et al. Replication Code.do
*	Date:  		July 19, 2018
*	Author: 	Michaël Aklin, S.P. Harish, Johannes Urpelainen
*	Data Used:  ReplicationData.dta
*	Output		na
*	Purpose:   	.do file to replicate the results of:
*				Aklin, Michaël, S.P. Harish, and Johannes Urpelainen. 
*				"A global analysis of progress in household electrification"
*				Energy Policy.
*
*				The code for the results listed in the appendix are
*				available from the authors.
*
*				To replicate our results, please modify the code below
*				where indicated. Note that to replicate the maps (Fig 1 and 2), 
*				you need to download additional files. The links are provided 
*				below. 
*	************************************************************************

clear all

*	HERE: comment out and replace working directory with the one in which
*	you downloaded the replication package.
// cd "..."
use "ReplicationData.dta"

*	************************************************************************
* 	A. Data: summary stats
*	************************************************************************

*	Summary statistics
quiet estpost sum elecrate_total elecrate_rural elecrate_urban ///
	delec_total_lowess delec_rural_lowess delec_urban_lowess	///
	wb_popdensityK ///
	pop_urban_rate rents_ff loghydropowercapitaK polity2 loggdpcap time ///
		if elecrate_total!=. | elecrate_rural!=. | elecrate_urban!=., d
esttab using "SummaryStats.tex", ///
	cells("mean(label(Mean) fmt(2)) p50(label(Median) fmt(1)) sd(label(S.D.) fmt(1)) min(label(Min.) fmt(1)) max(label(Max) fmt(1)) count(label(Obs.) fmt(0))") ///
	nonumber label replace noobs	///
	substitute(m$ m\\$)	
eststo clear


*	************************************************************************
* 	B. Analysis
*	************************************************************************

*	************************************************************************
* 	B.1 Table 1: Electrification rates
*	************************************************************************

*	Table 1, Panel A
eststo: xtreg ahu_total_w99 time, fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_total_w99 time if region=="East/Southeast Asia", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_total_w99 time if region=="Latin America & Caribbean", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_total_w99 time if region=="Middle East & North Africa", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_total_w99 time if region=="South Asia", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_total_w99 time if region=="Sub-Saharan Africa", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg wb_total_w99 time if ahu_total_sample==1, fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace

esttab using "Table1_A.tex"	///
	,  booktabs label replace ///
 	nodepvars se(2) b(2) noconstant ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	substitute(m$ m\\$)	///
	compress nogaps scalars("fixed Country FE" "N_g $\#$ Countries") r2(2) sfmt(%9.0f %9.0f %9.0f %9.0f) ///
	mtitles("All" "East Asia" "Latin America" "Middle East/N Africa" "South Asia" "SSA" "All")	///
	mgroups("New Database" "WB", pattern(1 0 0 0 0 0 1) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))	
eststo clear

*	Table 1, Panel B
eststo: xtreg ahu_rural_w99 time, fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_rural_w99 time if region=="East/Southeast Asia", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_rural_w99 time if region=="Latin America & Caribbean", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_rural_w99 time if region=="Middle East & North Africa", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_rural_w99 time if region=="South Asia", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_rural_w99 time if region=="Sub-Saharan Africa", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg wb_rural_w99 time if ahu_rural_sample==1, fe  cluster(ccode)
estadd local fixed "$\checkmark$" , replace

esttab using "Table_1b.tex"	///
	,  booktabs label replace ///
 	nodepvars se(2) b(2) noconstant ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	substitute(m$ m\\$)	///
	compress nogaps scalars("fixed Country FE" "N_g $\#$ Countries") r2(2) sfmt(%9.0f %9.0f %9.0f %9.0f) ///
	mtitles("All" "East Asia" "Latin America" "Middle East/N Africa" "South Asia" "SSA" "All")	///
	mgroups("New Database" "WB", pattern(1 0 0 0 0 0 1) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))	
eststo clear

*	Table 1, Panel C
eststo: xtreg ahu_urban_w99 time, fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_urban_w99 time if region=="East/Southeast Asia", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_urban_w99 time if region=="Latin America & Caribbean", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_urban_w99 time if region=="Middle East & North Africa", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_urban_w99 time if region=="South Asia", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg ahu_urban_w99 time if region=="Sub-Saharan Africa", fe cluster(ccode)
estadd local fixed "$\checkmark$" , replace
eststo: xtreg wb_urban_w99 time if ahu_urban_sample==1, fe  cluster(ccode)
estadd local fixed "$\checkmark$" , replace

esttab using "Table_1c.tex"	///
	,  booktabs label replace ///
 	nodepvars se(2) b(2) noconstant ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	substitute(m$ m\\$)	///
	compress nogaps scalars("fixed Country FE" "N_g $\#$ Countries") r2(2) sfmt(%9.0f %9.0f %9.0f %9.0f) ///
	mtitles("All" "East Asia" "Latin America" "Middle East/N Africa" "South Asia" "SSA" "All")	///
	mgroups("New Database" "WB", pattern(1 0 0 0 0 0 1) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))	
eststo clear

*	************************************************************************
* 	B.2. Table 2: Explaining discrepancies
*	Two outliers removed: Singapore and Timer-Leste
*	************************************************************************

*	Table 2, Panel A
eststo: reg delec_total_lowess wb_popdensityK time if outliers==0, cluster(ccode) 
eststo: reg delec_total_lowess pop_urban_rate time if outliers==0, cluster(ccode) 
eststo: reg delec_total_lowess loghydropowercapitaK time if outliers==0, cluster(ccode) 
eststo: reg delec_total_lowess rents_ff time if outliers==0, cluster(ccode) 
eststo: reg delec_total_lowess wb_natres_rents time if outliers==0, cluster(ccode) 
eststo: reg delec_total_lowess polity2 time if outliers==0, cluster(ccode) 
eststo: reg delec_total_lowess wb_popdensityK pop_urban_rate loghydropowercapitaK rents_ff wb_natres_rents polity2 time if outliers==0, cluster(ccode) 
estat vif

esttab using "Table_1a.tex"	///
	,  booktabs label replace ///
 	nodepvars se(2) b(2) ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	substitute(m$ m\\$)	///
	compress nogaps r2(2) sfmt(%9.0f %9.0f %9.0f %9.0f) ///
	nomtitles nonotes 
eststo clear

*	Table 2, Panel B
eststo: reg delec_rural_lowess wb_popdensityK time if outliers==0, cluster(ccode) 
eststo: reg delec_rural_lowess pop_urban_rate time if outliers==0, cluster(ccode) 
eststo: reg delec_rural_lowess loghydropowercapitaK time if outliers==0, cluster(ccode) 
eststo: reg delec_rural_lowess rents_ff time if outliers==0, cluster(ccode) 
eststo: reg delec_rural_lowess wb_natres_rents time if outliers==0, cluster(ccode) 
eststo: reg delec_rural_lowess polity2 time if outliers==0, cluster(ccode) 
eststo: reg delec_rural_lowess wb_popdensityK pop_urban_rate loghydropowercapitaK rents_ff wb_natres_rents polity2 time if outliers==0, cluster(ccode) 
estat vif

esttab using "Table_1b.tex"	///
	,  booktabs label replace ///
 	nodepvars se(2) b(2) ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	substitute(m$ m\\$)	///
	compress nogaps r2(2) sfmt(%9.0f %9.0f %9.0f %9.0f) ///
	nomtitles nonotes 
eststo clear

*	Table 2, Panel C
eststo: reg delec_urban_lowess wb_popdensityK time if outliers==0, cluster(ccode) 
eststo: reg delec_urban_lowess pop_urban_rate time if outliers==0, cluster(ccode) 
eststo: reg delec_urban_lowess loghydropowercapitaK time if outliers==0, cluster(ccode) 
eststo: reg delec_urban_lowess rents_ff time if outliers==0, cluster(ccode) 
eststo: reg delec_urban_lowess wb_natres_rents time if outliers==0, cluster(ccode) 
eststo: reg delec_urban_lowess polity2 time if outliers==0, cluster(ccode) 
eststo: reg delec_urban_lowess wb_popdensityK pop_urban_rate loghydropowercapitaK rents_ff wb_natres_rents polity2 time if outliers==0, cluster(ccode) 
estat vif

esttab using "Table_1c.tex"	///
	,  booktabs label replace ///
 	nodepvars se(2) b(2) ///
	star(* 0.10 ** 0.05 *** 0.01) ///
	substitute(m$ m\\$)	///
	compress nogaps r2(2) sfmt(%9.0f %9.0f %9.0f %9.0f) ///
	nomtitles nonotes 
eststo clear




*	************************************************************************
* 	C. Claims made in text
*	************************************************************************

* Countries with less than 30% electrification rates and income above $3k
tab countryname if elecrate_total<30 & ny_gdp_pcap_kd > 3000 & ny_gdp_pcap_kd != .

tab countryname if elecrate_rural<30 & ny_gdp_pcap_kd > 3000 & ny_gdp_pcap_kd != .


*	************************************************************************
* 	D. Appendix
*
*	Please contact the authors for the code used to generate the results
*	in the supplementary appendix.
*	************************************************************************


*	************************************************************************
* 	E. Graphs
*	************************************************************************


*	************************************************************************
*	Figure 1 and 2
*	Map: country trend: total and rural electrification, 5 year ma
*		Created with:
*		Code: http://www.stata.com/support/faqs/graphics/spmap-and-maps/
*		Maps: http://thematicmapping.org/downloads/world_borders.php 
*		The maps were then simplified to reduce their file size. To do
*		so: http://mapshaper.org/ (simplified to 2%).
*
*		To replicate the maps: you need to download the maps from the link above
*		and place them in the same folder as the replication package.
*	************************************************************************

preserve

keep if year == 1980 | year == 1985 | year == 1990 | year == 1995 | year == 2000 ///
	| year == 2005 | year == 2010 | year == 2015

keep countryname year ma_elecrate_total ma_elecrate_rural ma_elecrate_urban 

rename countryname NAME

* Code for map. HERE: ensure that the path to the maps is orrect. 
cd "./TM_WORLD_BORDERS-0.3/"

* Saving a temporary version of the dataset
save "temp.dta", replace

* Creating shape files
shp2dta using "tmw", database(usdb) coordinates(uscoord) genid(id) replace

* Modifying one of the files to merge with "temp.dta"
use usdb.dta, clear

replace NAME = "Bahamas, The" if NAME == "Bahamas"
replace NAME = "Cabo Verde" if NAME == "Cape Verde"
replace NAME = "Congo, Dem. Rep." if NAME == "Democratic Republic of the Congo"
replace NAME = "Congo, Rep." if NAME == "Congo"
replace NAME = "Myanmar" if NAME == "Burma"
replace NAME = "Gambia, The" if NAME == "Gambia"
replace NAME = "Hong Kong SAR, China" if NAME == "Hong Kong"
replace NAME = "Iran, Islamic Republic of" if NAME == "Iran (Islamic Republic of)"
replace NAME = "Korea, Dem. Rep." if NAME == "Korea, Democratic People's Republic of"
replace NAME = "Kyrgyz Republic" if NAME == "Kyrgyzstan"
replace NAME = "Libya" if NAME == "Libyan Arab Jamahiriya"
replace NAME = "Macao SAR, China" if NAME == "Macau"
replace NAME = "Russian Federation" if NAME == "Russia"
replace NAME = "Slovak Republic" if NAME == "Slovakia"
replace NAME = "Vietnam" if NAME == "Viet Nam"
replace NAME = "Tanzania" if NAME == "United Republic of Tanzania"
replace NAME = "Macedonia, FYR" if NAME == "The former Yugoslav Republic of Macedonia"
replace NAME = "Egypt, Arab Rep." if NAME == "Egypt"
replace NAME = "Iran, Islamic Rep." if NAME == "Iran, Islamic Republic of"
replace NAME = "Korea, Rep." if NAME == "Korea, Republic of"
replace NAME = "Lao PDR" if NAME == "Lao People's Democratic Republic"
replace NAME = "Micronesia, Fed. Sts." if NAME == "Micronesia, Federated States of"
replace NAME = "Faeroe Islands" if NAME == "Faroe Islands"
replace NAME = "Moldova" if NAME == "Republic of Moldova"
replace NAME = "St. Helena" if NAME == "Saint Helena"
replace NAME = "St. Kitts and Nevis" if NAME == "Saint Kitts and Nevis"
replace NAME = "St. Lucia" if NAME == "Saint Lucia"
replace NAME = "St. Vincent and the Grenadines" if NAME == "Saint Vincent and the Grenadines"
replace NAME = "Venezuela, RB" if NAME == "Venezuela"
replace NAME = "Yemen, Rep." if NAME == "Yemen"

*	Note: the 1's are small islands or tiny countries like Vatican.
merge 1:m NAME using temp.dta
drop if _merge == 2
drop _merge

*	Flag countries for which we have coded data
gen oecd = 0
replace oecd = 1 if NAME == "Australia"
replace oecd = 1 if NAME == "Austria"
replace oecd = 1 if NAME == "Belgium"
replace oecd = 1 if NAME == "Canada"
replace oecd = 1 if NAME == "Czech Republic"
replace oecd = 1 if NAME == "Denmark"
replace oecd = 1 if NAME == "Estonia"
replace oecd = 1 if NAME == "Finland"
replace oecd = 1 if NAME == "France"
replace oecd = 1 if NAME == "Germany"
replace oecd = 1 if NAME == "Greece"
replace oecd = 1 if NAME == "Hungary"
replace oecd = 1 if NAME == "Iceland"
replace oecd = 1 if NAME == "Ireland"
replace oecd = 1 if NAME == "Israel"
replace oecd = 1 if NAME == "Italy"
replace oecd = 1 if NAME == "Japan"
replace oecd = 1 if NAME == "Korea, Rep."
replace oecd = 1 if NAME == "Luxembourg"
replace oecd = 1 if NAME == "Netherlands"
replace oecd = 1 if NAME == "New Zealand"
replace oecd = 1 if NAME == "Norway"
replace oecd = 1 if NAME == "Poland"
replace oecd = 1 if NAME == "Portugal"
replace oecd = 1 if NAME == "Slovak Republic"
replace oecd = 1 if NAME == "Slovenia"
replace oecd = 1 if NAME == "Spain"
replace oecd = 1 if NAME == "Sweden"
replace oecd = 1 if NAME == "Switzerland"
replace oecd = 1 if NAME == "United Kingdom"
replace oecd = 1 if NAME == "United States"
replace oecd = 1 if NAME == "Romania"
replace oecd = 1 if NAME == "Bulgaria"
replace oecd = 1 if NAME == "Hungary"
replace oecd = 1 if NAME == "Slovak Republic"
replace oecd = 1 if NAME == "Czech Republic" 
replace oecd = 1 if NAME == "Poland"
replace oecd = 1 if NAME == "Moldova"
replace oecd = 1 if NAME == "Albania"
replace oecd = 1 if NAME == "Russian Federation"
replace oecd = 1 if NAME == "Ukraine"
replace oecd = 1 if NAME == "Belarus"
replace oecd = 1 if NAME == "Uzbekistan"
replace oecd = 1 if NAME == "Kazakhstan"
replace oecd = 1 if NAME == "Georgia"
replace oecd = 1 if NAME == "Azerbaijan"
replace oecd = 1 if NAME == "Lithuania"
replace oecd = 1 if NAME == "Moldova"
replace oecd = 1 if NAME == "Latvia"
replace oecd = 1 if NAME == "Kyrgyz Republic"
replace oecd = 1 if NAME == "Tajikistan"
replace oecd = 1 if NAME == "Armenia"
replace oecd = 1 if NAME == "Turkmenistan"
replace oecd = 1 if NAME == "Estonia"

drop if oecd==1


*	Total electrification in 1990
spmap ma_elecrate_total using uscoord if year == 1990 & id != 145   ///
		,  id(id) legtitle(varlab)	///
		fcolor(Blues2) ///
		clmethod(custom) clbreaks(0 10 20 30 40 50 60 70 80 90 100)	///
		legorder(hilo) legend( label(2 "0-10%")  ///
		label(3 "10-20%") label(4 "20-30%")	///
		label(5 "30-40%") label(6 "40-50%") label(7 "50-60%") label(8 "60-70%") ///
		label(9 "70-80%") label(10 "89-90%") label(11 "90-100%") label(12 "No Data"))	///
		ndfcolor(gray)
graph export "Figure_1a.pdf", replace 

*	Total electrification in 2010
spmap ma_elecrate_total using uscoord if year == 2010 & id != 145   ///
		,  id(id) legtitle(varlab)	///
		fcolor(Blues2) ///
		clmethod(custom) clbreaks(0 10 20 30 40 50 60 70 80 90 100)	///
		legorder(hilo) legend( label(2 "0-10%")  ///
		label(3 "10-20%") label(4 "20-30%")	///
		label(5 "30-40%") label(6 "40-50%") label(7 "50-60%") label(8 "60-70%") ///
		label(9 "70-80%") label(10 "89-90%") label(11 "90-100%") label(12 "No Data"))	///
		ndfcolor(gray)
graph export "Figure_1b.pdf", replace 


*	Rural ectrification in 1990
spmap ma_elecrate_rural using uscoord if year == 1990 & id != 145   ///
		,  id(id) legtitle(varlab)	///
		fcolor(Blues2) ///
		clmethod(custom) clbreaks(0 10 20 30 40 50 60 70 80 90 100)	///
		legorder(hilo) legend( label(2 "0-10%")  ///
		label(3 "10-20%") label(4 "20-30%")	///
		label(5 "30-40%") label(6 "40-50%") label(7 "50-60%") label(8 "60-70%") ///
		label(9 "70-80%") label(10 "89-90%") label(11 "90-100%") label(12 "No Data"))	///
		ndfcolor(gray)
graph export "Figure_2a.pdf", replace
 
*	Rural electrification in 2010
spmap ma_elecrate_rural using uscoord if year == 2010 & id != 145   ///
		,  id(id) legtitle(varlab)	///
		fcolor(Blues2) ///
		clmethod(custom) clbreaks(0 10 20 30 40 50 60 70 80 90 100)	///
		legorder(hilo) legend( label(2 "0-10%")  ///
		label(3 "10-20%") label(4 "20-30%")	///
		label(5 "30-40%") label(6 "40-50%") label(7 "50-60%") label(8 "60-70%") ///
		label(9 "70-80%") label(10 "89-90%") label(11 "90-100%") label(12 "No Data"))	///
		ndfcolor(gray)
graph export "Figure_2b.pdf", replace 

restore

*	************************************************************************
*	Figure 3.
*		The original Figure 3 was created in R using ggplot. Here is the code 
*		that generates a similar figure in Stata.  
*	************************************************************************

*	Figure 3, left panel
tw (scatter ahu_total_w99 loggdpcap) ///
	(lowess ahu_total_w99 loggdpcap)	///
		, 	///
		legend(off)	///
		title("Total Electrification (%)")
graph export "Figure_3a.pdf", replace 

*	Figure 3, center panel
tw (scatter ahu_urban_w99 loggdpcap) ///
	(lowess ahu_urban_w99 loggdpcap) 	///
		, 	///
		legend(off)	///
		title("Urban Electrification (%)")
graph export "Figure_3b.pdf", replace 

*	Figure 3, right panel
tw (scatter ahu_rural_w99 loggdpcap) ///
	(lowess ahu_rural_w99 loggdpcap) 	///
		, 	///
		legend(off)	///
		title("Rural Electrification (%)")
graph export "Figure_3c.pdf", replace 

	
	

